#include <stdio.h>
#include <math.h>

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
//#include <iostream>
//using namespace std;
#define index_v(i,j) (i+(j)*iter_max)
    //char dreorth[] = "dreorth";
extern "C" void orth_mgs(float *A, int iter_max, int n, char *dreorth, float *V, int tid) {
    if (strcmp(dreorth, "dreorth") != 0) {
        dreorth = "nodreorth";
    }
    // printf("%s\n", dreorth);
    //std::cout<< dreorth<< std::endl;
    
    float *v = (float*)malloc(iter_max * sizeof(float));
    double alpha, nv;
    
    for (int i = 0; i < iter_max; i++) {
        V[index_v(i,0)] = A[index_v(i, 0)];
    }
   /* 
    double norm_v = 0.0;
    for (int i = 0; i < iter_max; i++) {
        norm_v += v[i] * v[i];
    }
    norm_v = sqrt(norm_v);
    
    for (int i = 0; i < iter_max; i++) {
        V[index_v(i,0)] = v[i] / norm_v;
    }
    */
    for (int k = 1; k < n; k++) {
        for (int i = 0; i < iter_max; i++) {
            v[i] = A[index_v(i, k)];
        }
        /*
            if(tid==12){
            printf("\nv=\n"); for (int i = 0; i < iter_max; i++) printf("%f ",v[i]);
            }*/
            /*
            if(tid==12){
            printf("\nVm=\n"); for (int i = 0; i < iter_max; i++)   printf("%f ",V[index_v(i,j)]);
            }*/
           /* 
            if(tid==12){
                printf("j=%d,k=%d,alpha=%f \n",j,k,alpha);
                printf("\nv=\n"); for (int i = 0; i < iter_max; i++) printf("%f ",v[i]);
            }*/
        for (int j = 0; j < k ; j++) {
            alpha = 0.0;
            for (int i = 0; i < iter_max; i++) {
                alpha += v[i] * V[index_v(i,j)];
            }
            
            for (int i = 0; i < iter_max; i++) {
                v[i] -= alpha * V[index_v(i,j)];
            }
        }
      
        if (strcmp(dreorth, "dreorth") == 0) {
            for (int j = 0; j < k ; j++) {
                alpha = 0.0;
                for (int i = 0; i < iter_max; i++) {
                    alpha += v[i] * V[index_v(i,j)];
                }
         //   printf("j=%d,k=%d,alpha=%f \n",j,k,alpha);
                
                for (int i = 0; i < iter_max; i++) {
                    v[i] -= alpha * V[index_v(i,j)];
                }
            }
        }
        nv = 0.0;
        for (int i = 0; i < iter_max; i++) {
            nv += v[i] * v[i];
        }
        nv = sqrt(nv);
        //printf("nv=%.2e",nv);
        
        if (nv <= 1e-15) {
            printf("\n orth_mgs: Breakdown, step %d, val=%g \n\n", k+1, nv);
            return;
        }
        
        for (int i = 0; i < iter_max; i++) {
            v[i] /= nv;
            V[index_v(i, k)] = v[i];
        }
    }
    
    free(v);
}



extern "C" void orth_vec(float *P, float *rR, int iter_max, int kR, int tid) {
    float norm_rR = 0.0;
    
    for (int i = 0; i < iter_max; i++) {
        norm_rR += rR[i] * rR[i];
    }
    norm_rR = sqrt(norm_rR);
    
    for (int i = 0; i < iter_max; i++) {
        rR[i] = rR[i] / norm_rR;
    }
    /*
    if(tid==12){
        printf("w=\n");
        for (int i = 0; i < iter_max; i++) {
            printf("%.2e ",rR[i]);
        }}
    
       if(tid==12){ 
        printf("P=\n");
        for (int j = 0; j < kR; j++) {
            for (int i = 0; i < iter_max; i++) {
                printf("%f ",P[index_v(i,j)]);
            }
        printf("\n");
        }
       }*/
    for (int j = 0; j < kR; j++) {
        float alpha = 0.0;
        
        for (int i = 0; i < iter_max; i++) {
            alpha += P[index_v(i, j)] * rR[i];
        }
        //if(tid==12)printf("alpha=%.2e\n",alpha);
        
        for (int i = 0; i < iter_max; i++) {
            rR[i] -= alpha * P[index_v(i,j)];
        }
    }
    
    for (int j = 0; j < kR; j++) {
        float alpha = 0.0;
        
        for (int i = 0; i < iter_max; i++) {
            alpha += P[index_v(i,j)] * rR[i];
        }
       // if(tid==12)printf("alpha=%.2e\n",alpha);
        
        for (int i = 0; i < iter_max; i++) {
            rR[i] -= alpha * P[index_v(i,j)];
        }
    }
    
    float norm_w = 0.0;
    
    for (int i = 0; i < iter_max; i++) {
        norm_w += rR[i] * rR[i];
    }
    norm_w = sqrt(norm_w);
    
    for (int i = 0; i < iter_max; i++) {
        rR[i] /= norm_w;
    }
}


/*
int main() {
   
    float P[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
    float rR[] = {7.0, 8.0, 9.0};
    float w[3];
    
    orth_vec(P, rR, 3, 2, w);
    
    printf("w = [%f, %f, %f]\n", w[0], w[1], w[2]);


    float A[9] = {1, 2, 3, 4, 5, 1, 7, 8, 9};
    int iter_max = 3; // 行数
    int n = 3; // 列数
    char dreorth[] = "dreorth";
    float *V =new float[iter_max*n]; // 存储向量
    
    // 调用orth_mgs函数
    orth_mgs(A, iter_max, n, dreorth, V);
    
    // 打印结果
    printf("Orthogonal Matrix:\n");
    for (int i = 0; i < iter_max; i++) {
        for (int j = 0; j < n; j++) {
            printf("%.2f ", V[index_v(i,j)]);
        }
        printf("\n");
    }

    
    return 0;
}
*/
